%--------------------------------------------------------------------------
% computes the evolution of perturbations to the basic state by
% time-stepping the linearized equations. 
%
% here we assume that the interface is non-deformable.
%
%--------------------------------------------------------------------------


function lin_soln = stab(k, p, c)

global asy_base

asy_base = 1;

if nargin == 1
    p = params;
end


b_soln = base(p);

eta = linspace(0, 1, p.N)';

if nargin < 3
    
%     c = 0.01 * ones(size(eta));
%     c = c / norm(c, 'inf');
    
    % use eigenmode as initial condition
    [ev, k, M] = comp_eigs(k, p);
    [ef, ev] = eigs(M, 1, 'SM');
    c = -ef / norm(ef, 'inf');
end

opts = odeset('AbsTol',1e-10, 'RelTol',1e-10);
lin_soln = ode15s(@(t, c) pde(t, c, k, b_soln, p), [0, b_soln.x(end)], c, opts);

if (nargin == 1)
    subplot(2,1,1);
    surf(eta, lin_soln.x, lin_soln.y');
    contourf(eta, lin_soln.x, lin_soln.y', 100, 'linestyle','none');
    xlabel('$z / h(t)$','interpreter','latex','fontsize',12);
    ylabel('$t$','interpreter','latex','fontsize',12);
    title('contours of $\hat{c}(z,t)$','interpreter','latex','fontsize',12);
    colorbar
    subplot(2,1,2);
%     semilogx(lin_soln.x, lin_soln.y(p.N,:), 'k', b_soln.x, lin_soln.y(1,1) * test, 'r+');
    semilogx(lin_soln.x, lin_soln.y(p.N,:), 'k');
    xlabel('$t$','interpreter','latex','fontsize',12);
    ylabel('$\hat{c}(h(t),t)$','interpreter','latex','fontsize',12);
end


function dc = pde(t, c, k, b_soln, p)

global asy_base

if ~asy_base
    h_bar = deval(b_soln, t, p.N+1);
    c_bar = p.beta + deval(b_soln, t, 1:p.N);
    dcbar_dn = comp_base_deriv(c_bar, h_bar, p);
    
    eta = linspace(0, 1, p.N)';
    z = eta * h_bar;

else    
    h_bar = -(1 + lambertw(-p.beta / (p.beta - 1) * exp(((-p.beta + p.delta * t) / (p.beta - 1))))) * (p.beta - 1);
    eta = linspace(0, 1, p.N)';
    z = eta * h_bar;
    
    v0 = (1 - p.beta) * (1 - 1 ./ h_bar);
    c_bar = p.beta + v0 + p.delta * (-1/2) * (1 - p.beta) * (p.beta + v0) * z.^2 / h_bar^2 + p.delta * p.beta * (1-p.beta) / 6;
    dcbar_dn = -p.delta * (1 - p.beta) * (p.beta + v0) * eta;
end


% if (k >= 1 + 0 / sqrt(p.delta))
%     w = p.Ma / 2 * k * (h_bar - z) .* exp(k * (z - h_bar)) * c(p.N);
% else
    A = w_coeffs(c(p.N), h_bar, k, p);
    w = (A(1) * z + A(2)) .* cosh(k*z) + (A(3) * z + A(4)) .* sinh(k*z);
% end

d_hbar = -p.delta * c_bar(p.N);

dc_dn = [
    0;
    (c(3:end) - c(1:end-2)) / 2 / p.dz;
    -h_bar * p.delta * (1 - 2 * c_bar(p.N)) * c(p.N);
    ];


d2c_dn = [
    (2 * c(2) - 2 * c(1)) / p.dz^2;
    (c(3:end) - 2 * c(2:end-1) + c(1:end-2)) / p.dz^2;
    (2 * c(p.N-1) - 2 * c(p.N) + 2 * p.dz * dc_dn(p.N)) / p.dz^2
    ];

dc = eta * d_hbar .* dc_dn / h_bar - k^2 * c - w .* dcbar_dn / h_bar + d2c_dn / h_bar^2;